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1) Preconditioning matrices for Chebyshev derivative operators in several space dimension. 

We enclose a paper by E. Rothman explaining and summerising the problem. A joint paper 
with D. Funaro is being completed. 

2) The Jacobi matrix technique in computational fluid dynamics. 

We enclose the Ph.D. thesis of T. Sharp that was completed under this grant. 

3) Chebyshev techniques for periodic problems. 

We have developed Chebyshev techniques for periodic problems. These techniques give rise 
to differentiation matrices that have purely imaginary eigenvalues for odd number of derivatives 
and real and negative for even derivatives. We have demonstrated that for flows with boundary 
layer behaviour the periodic Chebyshev method is superior to Fourier methods. 

Wai-Sun Don, a graduate student in the division, is applying those ideas to simulate flows around 
a circular cylinder. Dr. D. Rudi of the Computational Method Branch provides supervision for 
this project. 
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ARSTRACT: 

The problem of preconditioning the matrices arising from pseudo-spectral Chebyshev 
approximations of first order operators is considered in both one and two dimensions. In 
one dimension a preconditioner represented by a full matirx which leads to preconditioned 
eigenvalues that are real, positive and lie between 1 and n!2, is already available. Since 
there are cases in which it is not computationally convenient to work with such a 
preconditioner, we study a large number of preconditioners which are "more sparse" (in 
particular three and four diagonal matrices). The eigenvalues of such preconditioned 
matrices are compared. In particular, the analysis is carried out for the quantity 

maxi Xil/minI Xil, where X^ are the preconditioned eigenvalues. 

We apply the results to the problem of finding the steady state solution to an 
equation of the type u^ = u^ + f, where Chebyshev collocation is used for the spatial 
variable and time discretization is performed by the Richardson method. 

In two dimensions different preconditioners are proposed for the matrix which arises 
from the pseudo-spectral discretization of the steady state problem in the square 
A = {(x,y,); - Kx^l, -l«y«l} 

U, . Uy = f 

U(x, y, 0) = Uq 

with boundary conditions at x = 1 and y = 1. Results are given for the CPU time and 
the number of iterations using a Richardson iteration method for the unpreconditioned 
and preconditioned cases. 
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1, INTRODUCTION 

To obtain the pseudo-spectral or collocation approximation let Pj^ be an interpolation 
operator. Let f(x) be a sufficiently smooth function defined in [-1,1] where f(x)=0 at the 
appropriate boundaries which yields a well-posed problem for (1.1). Then Pj^f is the 
interpolation of f at the collocation points Xj, i.e. 

P^fCXj) = f(Xj) and P^fcB^. j = 0 N 

To obtain a Chebyshev Gauss-Lobatto pseudo-spectal approximation in the interval 
[-1,1] we choose Xj = cos jn/N (j = 0,...,N), which when j * 0,N are the extrema of the 
N order Chebyshev polynomials Tj^j(x) = cos(N cos'^x). In order to construct the 
interpolant of f(x) at x, we define the polynomials 

(l-x2)T^(x)(-l)j*i 

gj(x) = 0 = 0,...,N) 

CjN^(x - Xj) 


1 (l«jSN-l). 

One can easily see that gj(x^) = 5j^. 

The N* degree interpolation polynomial Pj^f to f is given by 
N 

(1.8) Pf^f(x) = I f()^)gj(x) xeR 

j=o 

We must now be able to express derivatives of Pj^f in terms of f at the collocation 
points Xj. Differentiating (1.8) we obtain 

d"Pj^f(x) N jjn 

(1.9) = X f(5q) gj(x) 

dx" j=o dx" 


so that 


dx” 


= I f(Xj) (D„)^j 

j=0 


( 1 . 10 ) 
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where 

(1-11) (DnV = Sr- SkW 

^ X = Xj 

The pseudo-spectral Chebyshev derivative operator can be represented by the N x N 
matrix S„ = [s^j}, 
where 

Cj(-l)hk 

Sj^ = 

C^(Xj - x^) 

-^j 

Sji = , 

2(1 - Xj2) 


(k *i) 

2 N^+ 1 

% " ^ = ‘^NN 


In particular the Chebyshev pseudo-spectral approximation for u^ = u^, u(x,0) = Uq(x) 


N 


is given by Uj^ = I u(Xj,t)gj(x) and 
k=o 


6 u, 


N 


N 


0 t ~^ 2 !^u(Xj,t)Sj,j . 
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2. OUTLINE OF THE PROBLEM 
In the interval [-1,1] let 
2j-l 

= cos n 0 = 1.2.-.N) 

J 2N 

lies between Xj and Xj_^ 

“1 o I 2 o I o_ 

% %-l ^j-1 


O ^ O 1 

ICl Xq 


The modes Xj and have the following properties; 

Tn( 5P = 0 j = l....,N, and Tj^(xp = (-1)* i = 0 N 

Then consider the pseudo-spectral Chebyshev derivative operator with homogeneous 
boundary conditions at x = 1. This operator can be represented by the N x N matrix 

Sn = ■tSij}’ 

-2N2 + 1 
6 

Ci(-l)i 

I 

Cj(Xi- Xq) 

■"i 

. 2(1 - xf) 

The matrix Sj^ is full. The condition number C(Sf^) of is large. We have the 

following result which was obtained by Daniele Funaro. 

Lemma 1-1 : The condition number of Sj^ increases at least like N^. 

Proof : Let U • U denote the norm in £(IR^JR*^). Then the condition number 

of is given by 

c(Sn) = ISf^B . 


if i = j = N 
if i j 


if i = j = 1,...,N-1 


( 2 . 1 ) 
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It is known that ISjql^p(S,^)^CjN^, where p(Sj^) is the spectral radius of Sj^ and 
is a constant independent of N. On the other hand, we have 


( 2 . 2 ) 


IISj,j^8 = sup ISJ,jV1|rN. 

= 1 


Choose 


*<>0 


1 


1 


j = Furthermore, 


. Then li<PoH |jf= 1 and (S^Vq)] = ~ 

''N 


(2.3) 


"SnVoV= n I 

J=i -1 


^ wdx = c^ 


where C 2 does not depend on N. 

This implies Finally, using (2.1), we get C(Sj.^)>C 3 N^. Ihis 

proves the claim. 

Although the condition number is particularly meaningful for numerical applications, 
its determination is generally very difficult. Another quantity which is meaningful for 
practical application is o(M) = m^l Xjl/min! X-1 where M is an N x N matrix and X^ 
i = 1,...,N are its eigenvalues. It can be shown empirically that a(Sj^) behaves like N. 

We are interested in finding a "preconditioner" for In particular we are 

concerned with finding a matrix Rj,j such that the quantitity a(RJ^^Sj,j) is small. In 
general this does not imply that the corresponding condition number will be small (so the 
word "preconditioner" is not correct). 

In [DF ] Sjj is preconditioned by Rj^ = Z^Dj,^ where the N x N matrbc Dj,, = {dy} 
is defined by 

dji = X;) i = 1,...,N 

‘*u-l = l/(Xj_i- Xj) i = 2,...N 


otherwise 
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Hence D is the upwind finite differences matrix relative to the grid Xj. is 

the operator which maps the values of a polynomial in Pj^.j at the staggered grid points 
{^ 1 , into the values at the mesh points {Xj,...;Cjg}. 

Preconditioning by Zj^Dj^ results in preconditioned eigenvalues that are real, positive 
and lie between 1 and 7T/2. The ratio o(DJ^^Zj^^Sj,j) is bounded by nH (see [DF]). 
This is particularly interesting when the solution of the system 
(2.4) S„u + f = 0 

has to be found. If an iterative method is used, iterating = (ZjjDj^)'^Sj^ instead of 
Sj^ results in convergence in a few iterations. In the end of the computation the system 
(Zj,jDj^)u + f = 0 has to be solved. So we require that the matrix = Zj^D^j can be 
inverted easily. Although Z is a full matrix, it can be inverted very inexpensively in N 
log N operations. Thus the matrix Rj^ = Zj^Dj^ can be inverted very inexpensively. 
Nevertheless there are cases in which it is not computationally convenient to work with a 
full preconditioner. 

In particular when using an implicit method to find the solution at time T > 0 of 
the equation 


Uj = S^u + f, 

(2.5) u(x,0) = Uq(x), 

u(l,t) = 0 0«t«T. 

If the implicit Euler method is used, iterates of the matrix (I + AtSj^) are considered. A 
good preconditioner for this matrix turns out to be the matrix (I + AtZj^Dj^,). 
Unfortunately, due to the fact that Zj^ is a full matrix, (I + AtZj^Dj^) cannot be inverted 
inexpensively. In this work we shall present a large number of preconditioners which can 
be applied to the situations illustrated above. 
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3. ANALYSIS IN ONE DIMENSION 

In order to have the matrix I +M which can be easily inverted we substitute 

Zfj by some suitable matrix which has a simpler form and which has to be regarded as 
an approximation of the operator related to Zj^. 

The first idea is the following. Take Zj^ = {Zy} such that 


(3.1) 




= .5 


Z; 


ii+1 


= .5 


^NN 


1 . 


Then is a tridiagonal matrix so tha I + At Zj^Dj^ is also a tridiagonal matrbc. It 

can be shown that the eigenvalues of Mj^ = (ZjjDj^)'^Sj^ take the following form, 
k-l 

= k I Xj k = 1,...,N 

j=0 


Thus, the eigenvalues of the preconditioned matrix are real and positive. Hence we 
have o(Mf^) = N. The choice of Zj^ as in (3.1) corresponds to shifting the values from 
the staggered mesh to the initial mesh by averaging the two neighbour values. Instead of 
this we can choose Zn = {Zjj} corresponding to interpolation by first order polynomials. 




Tliis leads to the following definition of the matrix Zn = {Zy} 


(3.2) 


Z;: = 






Xj - 


i = 1,...,N-1 


i = 1.....N-1 


We found emperically, that the eigenvalues of the preconditioned matrix are still real and 
positive, but the quantity cf(M(,j) is now worse than that of the previous case. 

Another simple preconditioner which this time is not of the form is defined 
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by Rn 

= {rjj} where 



O 

II 

i = 1,...,N-1 


•"nN ~ ^(^N-1 ■ 


(3.3) 

= l/(^-l - 

i = 1,...,N-1 


- Vf) 

II 


%1,N ^ ■ %) 



In tliis case we still get real and positive preconditioned eigenvalues. Namely, they take 
the form: 


X., = 


k sin tt/N 


k = 




L J 

where is approximately equal to 2.46. We have a(Mj^) = N-1. This is the best we 
tested using tridiagonal preconditioners. Up to now the improvements are poor so we 
have to consider better approximations of the matrix Zj^. 

We consider an approximation of the operator related to Zj^ by interpolation with a 
polynomial of degree 2. One possible choice is the following: Z = {Zjj} where 


(3.4) 


^11 (^1 ~ ' ^ 2 ^ 

^12 “ (^1 * ^ lV (^2 ■ ^ 1 ^ 

Vi = (^i - 5i)(Xi - i = 2,...,N-1 

k = (^ - i = 2,..., N-1 

= (^i - - ?i)) i = 2,...,N-1 

^N,N-l ^ (% ■ ^N^/(^N-1 ■ 

^NN " (% ■ ' ^N-l) 


Now Z is a three-diagonal matrix, hence ZD is a four-diagonal matrix. The numerical 
experiments performed up to N = 32 give the following results. The eigenvalues X of 
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Mjj = are in general complex. The have po.sitive real part greater than .98. 

Most of them are concentrated at 1, Figure 3.1 shows the behavior of the \’s for N = 
12 and N = 20. 



Figure 3.1 - Location of the preconditioned eigenvalues in the complex plane. 

The corresponding o(Mj^) is represented in Table 3.1 for various N. This time 
a(Mj^) is bounded with respect to N. 


N 


8 

2.602 

16 

2.724 

24 

2.758 

32 

2.770 


Table 3.1 - Case of Four-diagonal preconditioner 
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other possibilities for the tridiagonal matrbc Z were tested. For example one 

possible choice, analogous to that of the matrix in (3.1), is to take = {Zy} such that 

(3.5) z„ = 3/4 

= 3/8 

Among all the experiments, the matrix proposed in (3.4) gives the best results. 
Five-diagonal preconditioners were tested, through interpolation by third degree polynomials. 
The results do not imporve those corresponding to Table 3.1. 

Now, if the system of linear equations (1.5) is solved, for example, by the implicit 
Richardson method, we are concerned with a(Mf^), where = (I + At Zj^Dfj)‘^(I + At Sv.^) 
and 7 .^ is the matrix in (3.4). The graph of o(Mj^) versus At is reported in Figure 2.2 
for some values of N. So, in addition to the fast convergence of the iterative schemiC, 
the preconditioning matrbc can be inverted efficiently. 






2 3 4 


5 


-At 


Figure 3.2 
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Now, if the steady state solution of problem (2.5) has to be found (which 
corresponds to the solution of problem (2.4), the explicit Richardson method can be used. 
This leads to the iterative scheme 

(3.6) u""^^ = (I + ^tSfj)u" , ncN 

If X denotes the general eigenvalue of Sj^, the scheme is stable provided we choose At 
such that 

r -2ReX ) 

(3.7) 0 < At < inf . 

X^ iXl2 J 

Within the stability region we have p(I + AtSj^) < 1 where p denotes the spectral radius. 

In order to speed up the convergence we experimentally find At* such that 
p(I + At*Sj^) attains its minimum inside the interval of stability. In general At* is not 
available. We use it here only to compare preconditioners. The same experiments are 

made for the preconditioned schemes, i.e.; 

(3.8) u""l = (I + AtZ^Df^)-l(I + AtSj^)u" 

(3.9) u"*l = (I + AtZ^Dj^)-l(I + AtS^)u" 

where Zn and Zj^ are respectively, the full matrix proposed in [DF] and the tridiagonal 
matrk given by (3.4). 

Both the schemes (3.8) and (3.9) converge to the same solution of (3.6). The 
results of these experiments are reported in Tables 3-2, 3-3, and 3-4 


N 

Maximum At for Stability 

At* 

p Corresponding to At* 

8 

-05461 

.02724 

.9817 

16 

.01836 

.009120 

.9919 

24 

.009232 

.005244 

.9939 

32 

.005138 1 

.002521 

.9966 


Table 3-2 Minimum spectral radius of the unpreconditioned amplification matrix 
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N 

Maximum At for Stability 

At* 

p Corresponding to At* 

8 

1.281 


.7810 

.2190 

16 

1.275 


.7787 

.2213 

24 

1.274 


.7783 

.2217 

32 

1.274 


.7782 

.2218 

Table 3-3 Minimum spectral radius 

of the preconditioned matrix : case of Zj,j. 

N 

1 Maximum At for Stability | 

At* 

p Corresponding to At* 

8 

.7685 


.5552 

.4448 

16 

.6381 


.3190 

.7112 

24 

.4258 


.2129 

.8465 

32 

.3210 


.1605 

.9043 


Table 3-4 Minimum spectral spectral radius of the preconditioned matrix : case of Z^, 

We also considered the second order Runge-Kutta scheme. This leads to the iterative 
scheme 

At^ 

(3.10) u""l = (I + AtSj^ + — S2)u" . neN 

In this case the stability restriction on At is given by 

(3.11) At^lXl'^ + 4 At^(ReX)|X|2 + 8At(ReX)2 + 8ReX < 0. 


for all eigenvalues, X, of We obtained results that were qualitatively analogous to 

that of Tables 3-2, 3-3, and 3-4. 
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4, ANALYSIS IN TWO DIMENSIONS 

In two dimensions we will consider the following steady state problem on the square 

A = {(x,y); -1 ^ X ^ 1, -1 $ y « 1} with homogeneous boundary conditions at x = 1, 

and y = 1. 

(4.1) u^ + Uy = f 

u(x,y,0) = Uq (x,y)cA 

u(l,y,t) = u(x,l,t) = 0. t ^ 0 

Tlie essential idea in obtaining a pseudo-spectral approximation to (4.1) is the same as it 
was in 1-dimension. That is, to approximate spacial derivatives by constructing a global 
interpolant through discrete points. To obtain the Chebyshev pseudo-spectral approximation 
we take as these points Xj = yj = cos rzj/N for j = 1,2,...,N. This means that we must 
interpolate at the points (Xj,yj) for i,j = 1,2,..., N. Consequently, the Chebyshev 

derivative operator for this problem can be represented by the x matrix + 

P'Sf^<2)p, where is a block diagonal matrix whose blocks are each equal to Sj,j, and P 
is a permutation matrix. If one orders the points (Xj,yj) by rows then 
corresponds to the derivative in the x-direction and P is constructed so that P^S^,^^P 
corresponds to the derivative in the y-direction. Without preconditioning + P^ S^^^P 
is ill-conditioned. 

As we saw in section 3 ZD is a good preconditioner for Sj,j. Thus, a natural 

approach to finding a preconditioner for + P^ S^^^P is to try + P'Z^^^D^^^P, 

where and are x block diagonal matrices whose blocks are the N x N 
matrices Z and D, respectively. 

To analyze the behavior of the eigenvalues of the preconditioning matrix, we define X 
as the generic eigenvalue of the preconditioned matrix, pj^ = m^l Xjl/minl \j|(i =1,2,...N^), 

is the maximum o such that ReX^o, and is the minimum r such that lx-ll«r. 
In particular, rj,j and give us an idea of the location of the eigenvalues. (See figure 
4-1 below.) 
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Numerical experiments performed for N = 4,6,8 are summarized in Table 4-1 


N 




4 

1.530 

.530 

1.000 

6 

1.552 

.589 

1.000 

8 

1.560 

.622 

1.000 


Table 4-1 


Although the eigenvalues of this preconditioned matrix, + 

P' S^2)p)^ are well behaved the matrix is full and thus difficult to invert. Another 
approach to constructing a preconditioner is to substitute in + P^ the 

tridiagonal matrix Z defined by (2.4) in place of Z. We will denote this new x 
matrix by + P‘ This matrix represents a finite difference scheme 

depending on seven points as illustrated by the stencil in Figure 4-2. 

• (i + Ij) 

• • • t 

(ij-2) (i,j-l) (i,j)(i,j + l) 

• (i-lj) 

• (i-2j) 


Figure 4-2 
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Results similar to those presented in Table 4-1 are prsented in Table 4-2 for 


N 

Pn 

■"N 


4 

2.206 

.980 

.897 

6 

5.026 

1.453 

.488 

8 

8.494 

1.602 

.306 


Table 4-2 

Although pj^ corresponding to increases more quickly than p^i 

corresponding to + P^Z^^^D^^V, the matrix can be inverted more efficiently. This 

is because + P‘ Z<^>0(2)p is a banded matrix (with N lower codiagonals and 2N 

upper codiagonals). 

Another preconditioner that we considered was of the form 
2(2)pt2;(2)p(D(2) + ptD<%). 


As in the 1-dimensional case, if the steady-state solution is to be found, the explicit 
Richardson method can be used. We experimentally find At* such that p(I + At*Wj^,) 
attains its minimum inside the region of stability. The same experiments are made for the 
preconditioned matrices. The results of the experiments are reported in Tables 3-3, 3-4, 
and 3-5. 


N 

Maximum At for Stability 

At* 

p Corresponding to At* 

4 

.1509563 

.07547812 

.9248009 

6 

.05228855 

.02614427 

.9715761 

8 

.0273053 

.01365251 

.9817108 

10 

.01944351 

.009721750 

.9810511 


Table 4-3 Minimum spectral radius of the unpreconditioned amplification matrix 
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N 

Maximum At for Stability 

At* 

p Corresponding to At* 

4 

.9879928 

.7172827 

.6693591 

6 

.9094622 

.6584506 

.8027519 

8 

.7200746 

.4954113 

.8715243 

0 

.6038381 

.3623029 

.8995660 


Table 4-4 Minimum spectral radius of the preconditioned matrix: 


Case of Z^P‘ Z^P(D^ + P' D^P) 


N 

Maxiumum At for Stability 

At* 

p Corresponding to At* 

4 

1.009613 

.6946138 

.3763489 

6 

.815170 

.6798409 

.6681257 

8 

.7685055 

.6870439 

37895085 


Table 4-5 Minimum spectral radius of the preconditioned matrix: 

Case of Zj^D^ + P‘ Zj^.Dj^P 

We also considered the second order Runge-Kutta scheme, and we obtained results 
that were qualitatively analogous to those of Tables 4-3, 4-4, and 4-5. 

We applied the Richardson schemes in the unpreconditioned version and in the 
preconditioned versions using the preconditioners Z.j^Dj^ + P* Zj^D^^P and 
Zj^P^Zj^P(Dj^ + P^Dj^P), to find the solution of the model problem 

Uj = Ujj + Uy - osin(a(x+l)) + osin(a(y+l) 
u(x,y,0) = sin(y-l)sin(x-l) 
u(-l,y,t) = u(x,-l,t) = 0 

of the type in (4.1). We used the optimal At* listed in tables 4-3, 4-4, and 4-5. Using 
the L 2 norm we considered the scheme to converge when the exact error stabilized. We 
also calculated global CPU times. The experiments were performed on the IBM 3081 and 
the results are reported in Tables 4-6,4-7, and 4-8. 
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N 

No. of Iteration 
for Convergence 

Error 

CPU Time for Convergence 
(in seconds) 

4 

122 

.12405 X lO'i 

.21 

6 

400 

.104182 X 10'3 

1.33 

8 

900 

.47112 X 10-^ 

5.82 


Table 4-6 Unpreconditional Euler 


N 

No. of Interation 
for Convergence 

Error 

CPU Time for 
Convergence 

CPU Time for Construction 
of Preconditioner and 
Finding LU Decomposition 
of Preconditioner 

4 

11 

.1246626 X 10-1 

-03 

.00 

6 

34 

.10400 X 10 -3 

.26 

.04 

8 

69 

.4709596 X 10- ® 

1.17 

.11 


Table 4-7 Preconditioned Richardson Case of Zj^D + P' 


N 

No. of Interations 
for Convergence 

Error 

CPU Time for 
Convergence 

CPU Time for 
Construction of 
Precondtioner and 
Finding LU 
Decomposition of 
Preconditioner 

4 

32 

.12468 X lO'l 

.09 

.00 

6 

173 

-104255 X 10'3 

1.10 

.00 

8 

263 

.4709613 X 10-6 

3.31 

.00 


Table 4-8 Precondition Richardson Case of Zj^P^ Zj^P(D + P* DP). 
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INTRODUCTION 


Flow past a body is, in general, specified by a variety of parameters such as 
thickness, angle of attack, camber, Mach number etc. A particular flow is, therefore, 
characterized by a single point in the corresponding parameter space. Conversely, the 
numerical calculation of a particular flow field yields information at just one point of 
the parameter space. However, the nature of a continuous range of nearby flow fields 
is of fundamental significance in the design and performance of aircraft. To treat this 
generally, one can consider the variational equations (which are linear) obtained by 
differentiating the exact equations with respect to each of the relevant parameters. The 
resulting matrix of derivatives of flow quantities is referred to as the Jacobi matrix. 

The subsequent procedure is, in principle, straightforward. One integrates the 
nonlinear governing equations - which results in the determination of just one point in 
parameter space - and simultaneously the variational equations governing the Jacobi 
matrix. The last is then used to describe the neighborhood of the already determined 
point of the parameter space. A method is presented herein which allows efficient 
generation of solutions in the neighborhood of a base solution. Since the variational 
equations are linear, the additional computational time required for their integration is 
modest. 

We have applied the Jacobi matrix technique to the direct calculation of 
inviscid supersonic flow about 

o two dimensional airfoils of varying thickness, angle of attack and camber 
o axisymmetric bodies of varying thickness and taper 
and the design (inverse) calculation of inviscid supersonic flow past 

o airfoils described by a given family of pressure distributions 
o axisymmetric bodies described by a given family of pressure distributions. 
Also we applied the method to subsonic potential flow about two dimensional 
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airfoils by modifying Jameson’s FL036. 

Results of our calculations show that the Jacobi method allows for the 
efficient and accurate generation of parametric solutions in the neighborhood of a 
known solution. In general terms, we consider a system of nonlinear partial differential 
equations 

E = 0 (LI) 

in the flow variables u, dependent variables x, and parameters g. For purposes of 
exposition we regard u(x;£) as known and seek the solution at a neighboring point in 
parameter space. The parametrically differentiated dependent variables are governed by 
the equations obtained by differentiating (1), viz. 


^ Fj(u,x;€)= 9E + M = 0 


( 1 . 2 ) 


The, in general, non-square matrix is known as the Jacobi matrix and the above 

36 ^ 

procedure provides a linear system of equations governing the Jacobi matrix. The term 
3Fj/3Uj actually represents an operator, the details of which are best left to the 
individual cases. If u°= u(x;ig) represents a known solution of the flow then any 
neighboring flow at some fixed point x is determined by 


o 3u° 

u(x;i)«*u°+ ^(i^e- Lo) 


(1.3) 


In what follows we will be somewhat loose in not distinguishing between the two sides 
of (1.3). A basic difficulty with what has been just said, in particular to the use of 
(1.3), is the fact that the conditions on the problem occur at locations which vary with 
€.. Specifically, both the boundary locations (and shock locations) may vary with changes 
in the parameters e.. We first present a method that avoids the difficulties implicit in 
such spatial variations with and later treat directly the formulation implicit in 
(1.1-3). 
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Chapter I 

The Jacobi Matrix Technique and 
It’s Application to Two-Dimensional 
Supersonic Flow 


The White Rabbit put on his spectacles. "Where shall I begin, please your 
Majesty?" he asked. 

"Begin at the beginning," the King said, very gravely, "and go on till you come to 
the end: then stop." 

- Alice’s Adventures in Wonderland 
Lewis Carroll 
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APPUCATTON TO 2D SUPERSONIC FLOW 


To illustrate this method we consider steady, inviscid, supersonic flows past two 
dimensional airfoils. For this purpose and in order to be specific, consider a family 
of profiles depending on three parameters (thickness, camber, and angle of attack). For 
completeness, we summarize the methods used in solving such flows [1], [2]. The 
equations are written in characteristic form as follows; 


S 3 = 0 


(1.1) 

CD 

+ 

II 

sin2M s- 

27 ^ 

(1.2) 

(9 - P(^t) )3 = 

(1 - tan 0 tanp) ~ ®a 

(1.3) 


Here the coordinates (a,/3) correspond to the streamlines, a = constant, and the 
characteristics, B = constant (Figure 1). 0 is the flow deflection angle, m is the Mach 
angle and s is the entropy. P(>r) is the Prandtl function given by 

P(ix)= X^tan'^ (x'^^tang) - g , X = ( 7 +l)/( 7 -l) . (1.4) 

An advantage to solving the above characteristic form of the equations is that it 
generates a body fit, shock fit coordinate system. We mention in passing that since the 
equations are exact, they are valid in the hypersonic flow regime so long as such real 
gas effects as disassociation and ionization can be ignored. 

The physical coordinates x,y satisfy the relations [1], [2] 

y „ = x^tan (0+4), y 3 = X 3 tan 0 ( 1 . 5 ) 

The transformation to (a, 6 ) coordinates leaves open two arbitrary functions and these 
are fixed so that the shock is along a = 3 and the airfoil is positioned along the line 
a = 0 (Figure 2). Appropriate boundary conditions at the body are 
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( 1 . 6 ) 


x((0,S) =B, y(0,B) = f(S,€). 6(0.8) = tan-l(f3(8, € )) . 

The Rankine-Hugoniot conditions govern the jumps in 6, 4 and s at the shock. 
Written in terms of the shock angle n, they are given by 

1 (M^- l)t an^n - 1 

tan6 = 

tanr? M^) tan^n (1.7) 


0 .2(l+^w)(l + iw) 

6 6 

sin 4 = 

(1+w) (1+0.2M ) - (l+-w) (1+ iw) 
6 6 


s = 2.5 In (l+-w) + 3.51n (1+ Iw) - 3.5 In (1+w) 
6 6 

where 


w = M^sin^h - 1 


The shock angle is related to the coordinates as follows 


tanh = dy / dx I = I 

shock ^8 shock 


( 1 . 8 ) 


(1.9) 


( 1 . 10 ) 


( 1 . 11 ) 


In the above we have assumed a perfect gas with constant specific heats and 
hence that 

p =p'^exp [(7- l)s] (1.12) 

It should be noted that this formulation eliminates the difficulty mentioned in the 
Introduction. Namely, by using the (oc,B) - coordinate system, a quantity such as 

5p 

^(«, 3 ;€) 

signifies the variation with €. at fixed a and 8. In particular it gives the variation of 
pressure say fixed at the body, a = 0, or at the shock, a = 8. This makes the 

integration of the differential equations significantly simpler. 


- 5 - 


VARIATIONAL EQUATIONS 


We are interested in solutions to these equations at points near a known solution. 
To pursue this we differentiate all of the above equations with respect to a typical 
parameter of interest. In keeping with the remarks at the close of the previous 
section, we emphasize that differentiation is with respect to with a and 6 held 
fixed. 


The mechanics of the differentiation are straightforward but tedious. We represent 
by capitalized variables the differentiated variables; 


30 

3s 

3x 

3v 

3g 


• S = 

• X = 

II 

II 


36 

3e 

3e 

de 


(1.13) 


when (1.4), (1.5), (1.6), and (1.8) are parametrically differentiated we obtain, 
Se = 0 


(1.14) 


_ s „cos2jr „ : _ o ,, 

®cx + 3 Y+P (g)M'^. s^= 0 

7 2 


(1.15) 


+ 


sec‘^0 tan 


e - (1 - tan 8 tan/r) ^3 0, 


a 


X o 0 ® ct ^3 

Pfi(#^) + tPp4(<^)#^3~ - " sec^4 tan0]T' + (1- tan0 tang) — (Xg X^j^) (1.16) 

X« ^a. 


tan( 0+ g) + Xjjj sec^(0+ g)(0 + Y) 


(1.17) 


Yg = Xg tan 0 + sec^0 (118) 

It should be noted that we have dropped the specification that e be a vector. This 

has been done for ease of exposition. This can be done without loss of generality. 
Variation with respect to each parameter can be treated separately, since only first 
order variations are being considered. 
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At the shock the parametrically differentiated equations are 


dn cos 


de (Xa+Xg)' 


—2 [y a e) ^ Y^(x „+x g)- X „(y„+y 3)' Xg +y g) ] 


(1.19) 


s = 


17.5 3.5 

+ 


3.5 


6 +7w 6 +w 1+w 


dw 

de 


( 1 . 20 ) 


0.2Afi + ^ wl - 0.2ru^wl fl+ iw 1 fo.2M^ - i - ^ wl 
1-3 18 J 'iirJ 


Y = 


3 18 '* dw 


A^ sin 2 /f 


— (1.21) 

de 


dw dn ->7 1 

where 57 = M^sin 2h ^ • A - I*0.2m2(1*w) - (Itlw) 


de L sin m 


d0 dn 
de 


(M^-l)tan^n -1 


(1+ M^)+ (1 + 2^ M ^) tan^n 

2 2 


( 1 . 22 ) 


[(1+2^±1 m^) + (1 +^ M^)tan^n] 
2 2 


1 [(M^- 1 ) sec n( 1+ ^ M 2) 


+ ( 1 + ^M^)tan^nJ - [ (M ^-l)tan^n- l] ( 1 + 2 ^ )sec ^n] | 


In the actual integration (1.21-1.25) are applied at the shock 
a = 6 

At the body a = 0 the appropriate equations are 

00 9fR(B. e) , 

X = 0, Y = f , 6 = cos^e 

0e 0€ 


(1.23) 


(1.24) 


In writing (1.24) we revert to the general case in which many parameters are 
being considered. At this point we can simultaneously numerically integrate the 
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non-linear system and the variational equations. The calculation of the base flow is 
second order accurate [2]. The calculation of the new flow is first order in space, 
second order in the parameters of interest. The calculation of the two flows is 
interleaved in that after the flow along 6 = constant is computed by the base code, 
the parametric code then calculates the exact derivatives in order to obtain the 
variational flow. 

RESULTS 

As we have already mentioned the method applies generally to many independently 
varying parameters. As a typical use of the variational quantities we use Taylor’s 
theorem to consider the change in pressure, 

Pnew ^ Phase '*■ (1-25) 


where Cj represents the various parameters with the differential coefficients calculated 
holding a and B fixed and the zero subscript denotes a reference or base calculation. 

For example, if just thickness is considered and denoted by say e, then at the 
body 

p p 1 o) “ p rr •"! 

new base '-Qe -*a=0,B ” base 9e ^x,y=f 


fdp fdf -> 

y^ftaT 


(1.26) 


'0y -'x, y= 

The second form exhibits the result obtained if variation in the physical plane is 
considered. 

In the numerical calculations that are discussed, we have taken for f a family of 
shapes given by 


y = 2€x(1-x) - X tan a + lOxe (x- 1) (x - ^ ) c 


(1.27) 


Thus, € is the thickness ratio based on chord, A the mean chord angle of attack, and 
c a scaling factor for the camber (shape) function. Figure 3 shows the effect of 
changing just the thickness (A,c = 0). Here we have plotted the pressure distribution 
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on the upper surface and the airfoil which, for aesthetic reasons, has the lower surface 
plotted as a reflection of the upper surface. Note that the method gives good 
agreement with the exact solution even when the new thickness is fifty percent greater 
than the base thickness. 

More generally we consider variations in all three parameters. Thus, the pressure 
relation at the body is 



(1.28) 


Figures 4 through 6 show the effect of changing various combinations of thickness, 
angle of attack, and camber. Here we see that, although the airfoil configurations are 
markedly different, there is very good agreement between the parametrically generated 
pressure distribution and the exact pressure distribution for the new airfoil. 


INVERSE CASE 

The method which has been presented also works as well on the inverse or design 
problem where the pressure on the body is known, but the shape of the body is to be 
determined. Using the Bernoulli equation and the perfect gas law one may show [1] 


H = sin 




7- 1 

expi—j (s + Inyp)) 


7-1 2 7-1 

1 + — M - exp( — (s + 1 n7p)) 

L 7 


■]'1 


(1.29) 


This when differentiated, yields 


- 


where 


( 7 - 1 )' 

27 


*y X 2 X 2 

(1^2^ ^ _ ( 7- 1 )'( 1 -^^'M ) exp(w) 


[i+ 2^M - exp(w)] sin^g 


[l+^M - exp(w)] sin^g 


d (tn p) 
d€ 


(1.30) 


7-1 

w = — y- (s+ln p) 
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Obviously at the airfoil we can no longer use (1.6) since we are hoping to 
determine the shape of the airfoil. Instead we must use (1.5) and (1.29). Therefore, 
the parametrically differentiated equations (1.24) must be replaced by (1.17), (1.18) and 
(1.30). The integration may now proceed as in the direct case [2]. The results of the 
variational calculations are presented in Figures 7 through 9. Notice that even for a 
20% change in the logarithm of the pressure (corresponding to a 50% increase in 
thickness), the difference between the exact airfoil shape and the computed shape is 
less than 4%. 
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HGURE CAPTIONS 


Fig. 1. 

Fig. 2. 
Fig. 3. 

Fig. 4. 

Fig. 5. 

Fig. 6. 

Fig. 7. 
Fig. 8. 
Fig. 9. 


Body, C+ characteristics, streamlines and C characteristics 
(dashed) in physical (x,y) plane, from [1]. 

Body, C+ characteristics and streamlines in (<x,B) plane. 

Pressure distribution on 10% and 15% thick airfoils at M = 2 and 

10% and 15% airfoils. 

Pressure distribution on 10% thick airfoil, uncambered at 0 degree 

angle of attack and cambered (c=0.1) at 5 degrees angle of attack 

at M = 2 along with respective bodies. 

Pressure distribution on 10% thick airfoil, uncambered at 0 degree 

angle of attack and canbered (c=0.2) at 5 degrees angle of attack 

at M = 2 along with respective bodies. 

Pressure distribution on 10% tliick airfoil, uncambered at 0 degree 

angle of attack and cambered (c=0.2) at 10 degrees angle of 
attack at M = 2 along with respective bodies. 

Inverse case: pressure distribution on 10% and 12% airfoils, M = 2, 

along with generated bodies. Dashed airfoil is computed shape. 

Inverse case: pressure distribution on 10% and 15% airfoils, M = 2, 

along with generated bodies. Dashed airfoil is computed shape. 

Inverse case: pressure distribution on 10% and 12% airfoils, M = 4, 

along with generated bodies. Dashed airfoil is computed shape. 
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Chapter II 

The Application of the Jacobi Matrk Technique 
to Axisymmetric Supersonic Flow 


"Curiouser and curiouser!" cried Alice (she was so much surprised, that for the 
moment she quite forgot how to speak good English). 

- Alice’s Adventures in Wonderland 
Lewis Carroll 
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AmJCATION TO AXISYMMKTRtC SUPERSONIC FLOW 


As another illustration of this method, we consider steady, inviscid, supersonic 
flows past axisymmetric bodies. For this purpose consider a family of profiles 
depending on two parameters, thickness and taper. As in Chapter 1, we shall 
summarize the methods used in solving such flows [1], [2]. The equations are written 
in characteristic form as follows: 


= 0 


tan 0 tanfi 


a 


(0 . P(^) ) = sjgM s- (a) - 

2y tan 9 + tan/r r 

fl ''6 

(9 - P(k) )fi = (1 - tan 0 tan/r) „ + tanM 

*ot r 


( 2 . 1 ) 

( 2 . 2 ) 

(2.3) 


Here the coordinates (a,0) correspond to the streamlines, oc = constant, and the C* 
characteristics, 6 = constant (Figure 1). 0 is the flow deflection angle, n is the Mach 
angle and s is the entropy. P(#i) is the Prandtl function given by 

P(n)= X^tan'^ (x'^tanji)- »r , X = (7 + 1)/ (y-l) . (2.4) 

As in the two-dimensional case, these equations are exact and are valid in the 
hypersonic flow regime so long as such real gas effects as disassociation and ionization 
can be ignored. The physical coordinates x,r satisfy the relations [1], [2] 

^ <x " r 3 = X 0 tan 9 (2.5) 

As with Chapter 1, the transformation to (oc,B) coordinates leaves open two 
arbitrary functions; these are fixed so that the shock is along a = B and the body is 
positioned along the line a = 0 (Figure 2). Hence, we have a body fit, shock fit 
coordinate system. The boundary conditions at the body are 

x(0,B) = B, r(0,B) = f(B,i), 0(P,B) = tan^(f 3 (B), €) . (2.6) 
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The Rankine-Hugoniot conditions govern the jumps in 0, il and s at the shock. 
Written in terms of the shock angle n, they are given by 


1 (M2-l)tan2r? - 1 

tan9 = 

tann + (1+^ M^) tan^n 


(2.7) 


0.2(l+^w)(l+lw) 

^ 6 6 
sin II = 

(l+w) (1+0.2M ) - (l+2w) (1+ i\v) 

6 6 

s = 2.5 In (l+2w) + 3.5ln (l+i\v) - 3.5 In (l+w) 
6 6 

where 

w = M^sin^i? - 1 

The shock angle n is related to the coordinates as follows 


( 2 . 8 ) 


(2.9) 


( 2 . 10 ) 


I 

Xoj + Xg shock 


tann = dr / dx I 


shock 


( 2 . 11 ) 


We have assumed a perfect gas with constant specific heats and hence that 


p =py exp [(7- l)s] 


( 2 . 12 ) 


It should be noted that this formulation eliminates the difficulty mentioned in the 
Introduction. For by using the (<x,$) - coordinate system, a quantity such as 

9p 

signifies the variation with e. at fixed a and B. In particular it gives the variation of 
pressure say fixed at the body or at the shock. This makes the integration of the 
differential equations significantly simpler. 
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VARIATIONAL EQUATIONS 


As in Chapter 1, we differentiate the governing equations with respect to the 
parameter of interest, keeping the coordinates a and 3 held fixed. The differentiation 
although straightforward is tedious. If we write 


00 0 s 

0 = > S = — 

0 e 0 € 


0 x 0 r 0 p 

X = — ’R = — (2.13) 

0 e 0 € 0 € 


and parametrically differentiate (2.1), (2.2), (2.3) and (2.5) we then obtain, 
S3 = 0 


e. 


7 r(tan9 + tan iiy 2 


r^c^Btan^jr 

+ 0 . -R = 0 

r(tan 0 + tanaf 

®3 f __§_®sec ^0 tan/rl ©- (1 - tan 0 tan/r) ^ 0 = 

Lx„ J “ 


r^tan 0 tan/r 


a 


a 


X|50_ . ^3 


P^(/r)Y3+ [Pj^^(g)n3- - P g . sec^n tanO + — sec^/r lY 


‘a 


X 3 


+ ( 1 - tan 0 tan/r) — X^^) + 


tanfi ^3^^"^ 


R 


R 


a 


(2.14) 


(2.15) 


(2.16) 


Rjjj = Xqj tan( 0 + g) + X„ sec^ (0 + ii)(e + Y) 


1 ^ = X 3 tan 0 + Xq sec^O 

At the shock the parametrically differentiated equations are 

*3)+ R3(x^ + X3) - Xj^(r„ + r3) 


dn 

de 


cos^ 


(2.17) 

(2.18) 


- X3(r«+r3) 


(2.19) 
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17.5 

3.5 

3.5 

dw 

6 +7w 

6+w 

1+w 

de 

0.2A[' 

U 2wl 

- 0.2fl+^ 

3 18 


( 


( 2 . 20 ) 


Y = 


sin2p[ 


3 18 J dw 

de 


( 2 . 21 ) 


where 


^ = M^sin 2rj ^ ^ = 1+0.2M2(1 +w) - (l+^w) (1+lw) 


d9 dn 
de 


- “^^4- A 

de L sin 4? 


(M^-l)tan^n -i 


(1+ M^)+ (l+^^M^)1an^n 


[(l+^M^) + (l+^M2)tan2)7] 


2 [(M^- 1) sec n( 1+ ^ M 2) 


( 2 . 22 ) 


+ ( 1 + ^ M ^ )tan^^] ■ [(M ^-l)tan^n-l] (1+ )sec ^n] | 


In the actual integration (2.18-2.22) are applied at the shock 

a = 3 (2.23) 

At the body a = 0 the appropriate equations are 

2L.= 0, R = f ^ , 0 = ^ cos^ 0 (2.24) 

3§_ 3€_ 

In writing (2.24) we revert to the general case in which many parameters are being 
considered. Now we simultaneously numerically integrate the non-linear system and the 
variational equations. The calculation of the base flow is second order accurate [2]. The 
calculation of the new flow is first order in space, second order in the parameters of 
interest. The calculation of the two flows is interleaved in that after the flow along 
3=0 constant is computed by the base code, the parametric code then calculates the 
exact derivatives in order to obtain the variational flow. 
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RESULTC 


In the numerical calculations discussed, we have taken for f a family of shapes 
given by 

r = 2€x(1-x) + 10x€ (X- 1) (X - i ) c (2.25) 

2 

Here we have taken e to be the thickness ratio based on chord, and c as a scaling 
factor for the taper function. 

Figure 3 shows the effect of changing just the thickness (c = 0). Here we have 
plotted the pressure distribution on the upper surface and the body which, for 
aesthetics, has the lower surface plotted as a reflection of the upper surface. Note that 
the method gives good agreement with the exact solution even when the new thickness 
is 50% more than the base thickness. 


Figure 4 shows the effect of changing a combination of thickness, and taper. Here 
we see that, although the body configurations are markedly different, there is very good 
agreement between the parametrically generated pressure distribution and the exact 
pressure distribution for the new body. 

INVERSE CASE 

The method which has been presented also works quite well in the inverse or 
design problem where the pressure on the body is known, but the shape of the body 
shape is to be determined. 

Using the Bernoulli equation and the perfect gas law one may show [1] 


fi = sin 




7 - 1 

exp(— j (s + Inyp)) 


7-1 2 
1+ — M 
2 


7-1 

(s + •nT'P)) 


l1 


(2.26) 


This when differentiated, yields 
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( 7 - 1 )^ 

27 


M^)exp(w) ( 7 - 1)^( 1+^ 


[i+ 2^M - exp(w)]' sin"/i 


S = 

, 2 ., 27 


2/ 1 ) exp(w) 


[l+^M - exp(w)] 
2 


sin^lJL 


Where 


d (In p) 
d€ 

(2.27) 


7-1 

w = -y (s+ln p) 


At the body, equation (2.6) is no longer valid since we are attempting to 
determine the shape of the body. Instead we must use (2.5) and (2.26). Therefore, 

the parametrically differentiated equations (2.24) must be replaced by (2.17), (2.18) and 
(2.27). The integration may now proceed as in the direct case [2]. The results of the 
variational calculations are presented in Figures 5 and 6. Notice that even for a 10% 
change in the logarithm of the pressure (corresponding to a 20% increase in thickness), 
tihe difference between the exact body shape and the computed shape is less than 1%. 
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FIGURE CAPTIONS 


Fig. 1 

Fig. 2. 
Fig. 3. 

Fig. 4. 

Fig. 5. 

Fig. 6. 


Body, C+ characteristics, streamlines and C characteristics (dashed) in 
physical (x,y) plane, from [1]. 

Body, C+ characteristics and streamlines in (a,6) plane. 

Pressure distribution on 25 % and 30% thick bodies at M = 6 and 
the respective bodies. 

Pressure distribution on untapered, 25% thick body and 0.10 taper, 
30% thick bodies at M = 6 and the respective bodies. 

Inverse case: Pressure distribution on 25% and 30% thick bodies, 

M = 4 along with generated bodies. Dashed body is computed shape. 

Inverse case: Pressure distribution on 25% and 30% thick bodies, 

M = 6 along with generated bodies. Dashed body is computed shape. 
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Figure 1 


COMPUTATIONAL PLANE 
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Chapter III 

The Jacobi Matrix Method for General Flows 


Here one of the guinea-pigs 
officers of the court. 


cheered, and was immediately suppressed 


by the 


- Alrce’s Adventures in Wonderland 
Lewis Carroll 
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DIRRCr DIFFERENCING 


The procedure outlined in Chapters 1 and 2 holds in much greater generality than 
we have considered. The Jacobi matrix technique could also be applied to unsteady 

flows and to viscous flows in three dimensions. However, the method as presented so 
far, has one possible drawback which was alluded to earlier - to obtain the Jacobi 

matrix we must analytically differentiate the relevant equations and boundary conditions. 
In this chapter we propose a procedure which will allow for the calculation of the 
Jacobi matrix by the use of differential approximations. The goal is to obtain the 
Jacobi matrix, and hence be able to calculate a range of solutions in parameter space, 
using the results obtained from solving the nonlinear system (1.1) at only two distinct 
values of £. This differential approach will be applied to the case of two dimensional 
supersonic flow considered in Chapter 1 and to two dimensional subsonic potential flow. 

In the Introduction we said that if u® = u(x ; £q) represented a known solution 
of the base flow then any neighboring flow at some fixed point x is approximated by 

u(x;e) u° + (€^ - €^) (3.1) 


9uS 

The obvious first order approximation for is 

% 



As 

^ fixed 


u(x;i.) - u® 

- id 


^ fixed 




0U® 

What this says is that to compute we can take the value of u at the location x 

in the base flow and subtract it from the value of u found from the perturbed flow 
(£ = + A£) at the same location, hi practice, this may require interpolation on one 

computational grid. 

This approach requires special attention at a boundary. In our approach both 
material boundaries and possible shocks are taken to be boundaries and both give rise 
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to locations which change with This would certainly be the case if we chose to 
vary the parameters of a body. 

To be more specific, we would like to be able to use the calculation of pressure 
in the base flow in order to compute the pressure on the new body. Thus the 
formulas (1.25), (1.26) are no longer applicable since they apply at a fixed field point. 
Therefore, to correct (1.26) we must include changes in location of the body due to 
changes in €,. In the interests of simplicity we specify a three dimensional body by 

y = f(x,z;e) (3.3) 


A typical quantity, say pressure, at the new body, which we will specify by is 
related to the old body ^ in the following way 

9p(xo;i<)) ap(xo;io) af 

p(2^;i) ^ P(Xo;iq) + — (3.4) 
where and 

h = (3-5) 

and 

2^ = (Xo,f(Xo,Zo;i),Zo) (3.6) 


Compare (3.4) with equation (1.26). 

Note that we have related 2^ to ^ by placing 2^ directly above ^ in the x-z 

plane. Other choices are possible and may be more appropriate in certain cases. 

Equation (3.4) in fact gives us the ability to compute the pressure at the new 

9p 

body, but requires knowledge of the differential coefficients . They can be 

obtained from 


9p(xo;io) 


p(2^oiiJ - Pfeo’’^) ■ 


ap(Xo;€o) 


9% 




9yo 


9lk 






(3.7) 
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to first order, or 


3€q 

9p0^,§^) p(2^;(e + € 0 )*^) ■ P(^‘> 

^ 

(3.8) 

It should be noted that the differential determination of the differential coefficient 

9p 9p 

requires not only calculation of the different flow fields, but also of . 

9e. ay^ 

9p 

Therefore, in the numerical calculation it is necessary to compute — — at the body. 

9yo 

This we do by interpolation. 

To illustrate these remarks we return to the case treated in Chapter 1. We 

consider two dimensional supersonic flow at thicknesses of 10% and 10.1% to calculate 
dx dy d0 dfi ds 

, , , , . Using equation (3.8) the resulting derivatives were 

de de de de de 

then used to compute the pressue distribution on a 15% thick airfoil (Figure 1). Note 
that this pressure distributuion compares favorably with that computed from using the 
Jacobi matrix generated by solving the differential equations (Chapter 1, Figure 3). The 
error between the two computations is less then 1%. 

2D SURSONIC FLOW 

As a second illustration we apply the Jacobi matrix technique to the potential 
equation for two dimensional compressible flow. The potential equation is derived by 
assuming inviscid, irrotational flow and is valid for subsonic flows and for low transonic 
flows when boundary layer effects can be neglected. 

Since we have implemented the Jacobi technique by modifying Jameson’s 
computer code FL036 we will summarize the derivation of the relevant equations and 
their solution [1],[2]. 

Under the assumption of irrotational flow we may introduce a velocity potential 
(p such that 


£ 

T 


) - 


3p(Xq; 


3i^ 


) af 


9yo 
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( 3 . 9 ) 


U = 0^ V = 0y 

The potential satisfies the quasilinear equation 


(a^ - u^) <t>^ - 2uv + (a^ - v^) = 0 (3.10) 


where a is the local speed of sound. When given the ratio of specific heats y, the 
stagnation speed of sound a^ and the local speed q = ^ ^2 the speed of sound is 

determined by 


2 y-^ 

^ - -T- 


(3.11) 


We will consider (3.10) for subsonic flows. (But see Figure 9 for a transonic 
case). 

At the body the flow must satisfy the tangentcy conditionn 


dp 

8n 


= 0 


(3.12) 


where n is the normal derivative and the Kutta condition - that the tangential velocity 
is bounded at the trailing edge. In the far field the potential approaches the potential 
of a vortex in compressible flow and a uniform stream. The density and pressure are 
determined by relations 

p7-i = Mia^ (3.13) 

and 

py 

P = — (3.14) 

7MI 


The coordinate system used for computation is generated by conformally mapping 
the exterior of the airfoil to the interior of the unit circle. The airfoil itself becomes 
the coordinate line r = 1 (Figure 2). 

Since the far field boundaiy condition must now be applied at r = 0, where the 
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potential becomes infinite, a reduced potential which removes this singularity is 
introduced by 

cos (9 + a) 

G = 4>- ^ ^ + E(0 + a) (3.15) 

Here a is the angle of attack and 27ZE is the circulation. 

If the modulus of the transformation from the physical plane to the circle plane 
is denoted by H then (3.10) becomes 

g 

(a2 - u2)Gee - 2uvrG,e + (a^ - y2)r — (rG) 

- 2uv(Gq - E) + (u^ - v2)rGj. + (u^ + v^) ( Hg + vH^) = 0 (3.16) 

The u and v are the velocity components in the 9 and r directions, respectively and 
are given by 

r(G 0 -E)- sin(9 + a) r^Gj.- cos(9 + a) 

u = , V = (3.17) 

H H 

The Neumann boundary condition (3.12) becomes 

G = cos(0 + a) at r = 1 (3.18) 

while the far field condition is 

G = E{9 + a - tan"^ [Vi . tan(0 + a) ] } at r = 0 (3.19) 

The circulation is determined by the Kutta condition which requires that the velocity be 
finite at the trailing edge of the airfoil. Here we have H = 0 and i>Q = 0 so (3.15) 
reduces to 

E = Gq - sina at r = 1, 0 = 0 (3.20) 

The details of the calculation of H and of the multigrid solution of (3.16)-3.20) 
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are not essential for our purposes and are discussed in references [1], [3], [4], [5]. The 
important point is that the transformation to the circle plane is conformal so that every 
airfoil in the physical plane is mapped to a circle and every physical flow is mapped to 
the interior of the circle. 

CALCULATION OF TEiE .TACOBI MATRIX 

The equations for computing the Jacobi matrix by finite differences were given by 

equations (3.4)-(3.8). In the subsonic case the only parameter changed was the 

airfoil thickness based on chord. Due to the construction of equation (3.16) the 

quantities which are of interest are the reduced potential G, and the metric H. To 

define the locations x in equation (3.2) we note that in the circle plane the points are 

spaced angularly (9) as 27T/(the number of grid points about airfoil) and radially (r) as 

l/(the number of grid points fi’om airfoil to far field). Therefore, it is natural to 

define the location x by the intersection of these lines. 

The variational flow was computed using essentially the same procedure which was 

used to calculate the pressure in the two dimensional supersonic flow case. In the 

transformed plane we first compute the flow about an airfoil of thickness and save 

the converged values of G and H. Next we compute the flow about an airfoil of 

thickness Cq. We use these computed values of G and H along with those from the 

dG dH 

run at tnickness 6, to compute and using (3.7). G and H for the 

^ de de 

variational flow at thickness € is computed using equation (3.4). 

RESULTS 

Figure 3 shows the results of a parametric calculation using a base airfoil of 
10% thickness based on chord and a second airfoil of 10.1% thickness to predict the 
pressure distribution on a 14% thick airfoil. It should be noted that there is very 
close agreement between the parametric calculation and the solution given by FL036. 

Figure 4 uses a 10% thick airfoil and 10.1% thick airfoil to calculate the flow 
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over 15% thick profile - that is, a profile which is 50% more than the base airfoil. 
Again the agreement is quite good. Figures 5 through 8 show the same calculations 
for flows at different Mach numbers. All show close agreement between the 
parametrically generated solutions and those given by FL036. 

The method breaks down when there is a drastic change in the behavior of the 
solution in the parameter space. This is illustrated in Figure 9. Here the flows about 
the 10% and 10.1% thick airfoils are subsonic but the flow about the 15% thick airfoil 
is supercritical. The method is unable to account for the shock. 
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FIGURE C APTIONS 


Fig. 1. 

Fig. 2. 
Fig. 3. 

Fig. 4. 

Fig. 5. 

Fig. 6. 

Fig. 7. 

Fig. 8. 

Fig. 9. 


Pressure distribution on 10% and 15% thick airfoils at M = 4 
calculated by direct differencing along with respective bodies. 

Computational plane, from [2]. 

Pressure distribution on 14% thick airfoil, M = 0.75. Solid curve is 
FL036 result, dashed is parametric. Base and new airfoils are also 
shown. 

Pressure distribution on 15% thick airfoil, M = 0.75. Solid curve is 
FL036 result, dashed is parametric. Base and new airfoils are also 

shown. 

Pressure distribution on 14% thick airfoil, M = 0.60. Solid curve is 

FL036 result, dashed is parametric. Base and new airfoils are also 

shown. 

Pressure distribution on 15% thich airfoil, M = 0.60. Solid curve if 

FL036 result, dashed is parametric. Base and new airfoils are also 

shown. 

Pressure distribution on 14% thick airfoil, M = 0.45. Solid curve is 

FL036 result, dashed is parametric. Base and new airfoils are also 

shown. 

Pressure distribution on 15% thick airfoil, M = 0.45. Solic curve is 

FL036 result, dashed is parametric. Base and new airfoils are also 

shown. 

Pressure distribution on 15% thick airfoil, M = 0.80. Solid curve is 

FL036 results, dashed is parametric. Base and new airfoils are also 
shown. 
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AIRFOIL (r=1) 



F igure 2 
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